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Abstract Electrophoretic separation of a mixture of chemical species is a fun- 
damental technique of great usefulness in biology, health care and forensics. In 
capillary electrophoresis the sample migrates in a microcapillary in the presence 
of a background electrolyte. When the ionic concentration of the sample is suffi- 
ciently high, the signal is known to exhibit features reminiscent of nonlinear waves 
including sharp concentration 'shocks'. In this paper we consider a simplified model 
consisting of a single sample ion and a background electrolyte consisting of a sin- 
gle co-ion and a counterion in the absence of any processes that might change 
the ionization states of the constituents. If the ionic diffusivities are assumed to 
be the same for all constituents the concentration of sample ion is shown to obey 

_C ' a one dimensional advection diffusion equation with a concentration dependent 

advection velocity. If the analyte concentration is sufficiently low in a suitable 
non-dimensional sense, Burgers' equation is recovered, and thus, the time depen- 
dent problem is exactly solvable with arbitrary initial conditions. In the case of 
small diffusivity either a leading edge or trailing edge shock is formed depend- 
ing on the electrophoretic mobility of the sample ion relative to the background 

\f-*. . ions. Analytical formulas are presented for the shape, width and migration veloc- 

ity of the sample peak and it is shown that axial dispersion at long times may be 
characterized by an effective diffusivity that is exactly calculated. These results 

r ^J ■ are consistent with known observations from physical and numerical simulation 

f*^ , experiments. 
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1 Introduction 

Slab gel electrophoresis has long been a familiar technique in molecular biology for 
sorting biomolecules according to their size and charge. For example, the Sanger 
method of sequencing DNA and RNA uses gel electrophoresis for sorting the prod- 
ucts of the chain termination reaction. The SDS-PAGE technique used for sorting 
proteins by size uses electrophoresis in a polyacrylamide gel to separate the surfac- 
tant coated amino acid chains resulting from the denaturing of the proteins. The 
technique of Capillary Electrophoresis (CE) was introduced about thirty years 
ago and has largely replaced gel electrophoresis in many modern applications. In 
CE, the sample migrates through a microcapillary (about 50 /jm in diameter and 
tens of centimeters long) instead of a slab of porous gel. The micro-capillary is 
connected at the inlet and at the outlet to reservoirs, as shown in the schematic 
sketch in Figure [1] and the entire system is filled with an electrolytic buffer [1] . A 
large electrical voltage (tens of kilovolts) is applied with electrodes between the 
inlet and outlet reservoirs. The sample separates into distinct zones based on the 
migration speeds of individual components in the applied electric field. A detector 
placed near the outlet records a series of peaks (an electropherogram) that may be 
used for analyzing the sample. Often (but not always) there is an electroosmotic 
flow 2 present in the capillary that causes both positive and negative ions to be 
swept past a single detector placed near the outlet. For linear polyelectrolytes (e.g. 
DNA) the mobility is insensitive to the polymer length. In such cases a sieving 
medium such as a polyacrylamide gel may be loaded into the micro-capillary. 

CE has many advantages over the traditional slab gel electrophoresis method, 
one of the most important ones being ease of parallelization. Since many CE chan- 
nels can be etched on a single substrate, multiple separations can be run simul- 
taneously. Such scalability is critical to applications such as genome sequencing, 
and indeed CE is the method of choice in most such applications. Though CE is 
generally considered superior to the more traditional gel electrophoresis methods 
one of its limitations is the high demands placed on the sensitivity of the photo 
detector. The detector, which consists of a UV or laser source coupled to a photo 
cell must be sensitive to light attenuation over very short optical path lengths (the 
capillary diameter). Therefore, in order to maximize the likelihood of detection, 
it is preferable to use high sample concentrations. This requirement however con- 
flicts with the other requirement of CE which is to minimize peak broadening or 
dispersion in order to improve resolution and signal strength |3]. The conflict is 
caused by a phenomenon known as "electromigration dispersion" or the "sample 
overloading effect" . The underlying physical mechanism is well understood and 
may be roughly explained as follows: when the concentration of sample ions is 
non-negligible compared to that of the background electrolytes the local electrical 
conductivity (a) is altered in the vicinity of the peak. On the other hand, since the 
current (J) is constant, the electric field (E) must change locally since by Ohm's 
law J = oE. This varying electric field alters the effective migration speed of the 
sample (or analyte) ion which in turn alters its concentration distribution thereby 
giving rise to a nonlinear transport problem that must be solved in a self consis- 
tent manner. Electromigration dispersion causes highly asymmetric concentration 
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profiles, high rates of dispersion and shock like structures that are reminiscent 
of nonlinear waves seen in many other physical systems. These effects have been 
widely reported in the literature on electrophoresis (see e.g. @]). 

Simple one dimensional mathematical models of electromigration dispersion 
have been constructed 15,6,7,8 by invoking the assumption of vanishing diffusiv- 
ity of all species and the requirement of local charge neutrality, which for a three 
ion system yields a single nonlinear hyperbolic equation. Solutions describe the 
observed steepening of an intitially smooth profile leading to subsequent shock 
formation. Two recent reviews [9lll0j provide more extensive references to one 
dimensional models and the behavior of their solutions as evident from numeri- 
cal simulations. These models assume that only strong electrolytes are involved 
so that the valence states of the constituent ions are fixed. This assumption is 
inadequate when ampholytes (such as proteins) are present. The separation of 
ampholytes is usually performed by exploiting two other modes of electrophoresis: 
Isotachophoresis and Isoelectric focusing. The mode of electrophoresis described 
in the introduction and discussed in this paper by contrast is called "Zone elec- 
trophoresis" . From a mathematical perspective the various modes of electrophore- 
sis correspond to different kinds of initial and boundary conditions of the same 
governing equations. An excellent discussion of the mathematical problem involved 
in the various modes of electrophoresis may be found in Chapter 3 of the book by 
Fife [11] and the practical aspects are discussed in a number of text books [III 12] . 
The essential technique for reducing the complexity of the transport equations in 
the presence of ampholytes was developed by Saville and Palusinski [13_ using the 
well founded assumption that the time scale associated with ionization phenomena 
is very short compared to the time scale of diffusive or electrophoretic transport. 
Thus, local equilibrium may be assumed to exist between the ampholyte and its 
dissociation products so that the concentrations of the various charged states can 
be algebraically related to the concentration of the ampholyte in the neutral state. 
A historically contentious issue and an area of some mathematical interest concerns 
the assumption of local charge neutrality mentioned earlier. Local charge neutrality 
results when the small parameter (in some appropriate dimensionless formulation) 
multiplying the Laplacian of the potential in Poisson's equation vanishes. This as- 
sumption, first invoked by Planck 14 , results in a singular perturbation problem 
leading to boundary layers and associated interesting physics in various areas in- 
cluding semiconductor junctions and in selectively permeable membranes [15] - In 
the context of electrophoresis, the singular nature of the local charge neutrality 
assumption is manifested at the sharp boundaries of isotachophoretic shocks where 
the shock is regularized by finiteness of the small parameter. This aspect has been 
studied by Fife [16] who also presented some results relating to the regularity of 
solutions. A proof of the global existence of solutions for the coupled species and 
electric field equations in one dimension in the absence of ampholytes have been 
presented by Avrin [T7]. In the present paper, we consider the one dimensional 
model of electromigration dispersion in the absence of ampholytes, however, we 
replace the assumption of zero ionic diffusivities with one of equal diffusivity of 
the ionic constituents. We show that this simple modification leads to a nonlinear 
advection-diffusion equation which reduces to Burgers' equation for the analyte 
concentration if this concentration is not too large. The wealth of information that 
is available about the Burgers' equation can thereby be exploited for understanding 
electromigration dispersion. 
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2 Fundamental equations of transport 

We consider a system of N ionic species and denote by Ci(x,t) the number of 
ions of species i (i = 1, 2, . . . , N) present per unit volume. Here x is the direction 
along the capillary and t is time. The system is presumed one dimensional, that 
is, the dependent variables are assumed not to vary over the cross-section of the 
capillary. In particular, this implies, electroosmotic flow, if present, is of the "plug 
flow" type. Cross stream variations can and do arise in CE and have been widely 
studied but it is not central to the present discussion. The concentrations obey 
the equations 

§ + §! = » w 

where qi is the flux (amount passing per unit area of cross section per unit time) 
of species i. In writing Eq. ([l|) we have presumed that there are no sources or sinks 
for the ions. In particular, ionization and recombination phenomena are ignored, 
which is justified if only strong electrolytes are involved. The fluxes q^ are related 
to the concentration fields through the Nerst-Planck model 

q i (x,t) = -D i -^ + n i c i (x,t)E(x,t), (2) 

where E(x,t) is the electric field (directed along the capillary), m is the elec- 
trophoretic mobility and Di is the diffusivity of the 'ith' species. In addition to 
the electrophoretic mobility, which is the velocity acquired by the ion per unit of 
applied electric field, we would also refer to the mobility (ui) which is the veloc- 
ity per unit applied force. The two are relateqj to each other as Hi = zieui and 
related to the diffusivity through the Stokes-Einstein relation [18] Dj/ttj = kgT, 
where kg is Boltzmann's constant, T is the absolute temperature, e denotes the 
proton charge and z.j is the signed valence of the ith species which can be posi- 
tive, negative or zero. When fluid flow with velocity V is present, a term Vci is 
added to the right hand side of Eq. ([2}. However, if the flow is an electroosmotic 
plug flow, this term may be absorbed in the second term by simply redefining ^ 
as the "total" instead of the "electrophoretic" mobility. If the capillary is filled 
with a polymer network then the mobility and diffusivity introduced above should 
be interpreted as the "effective" values of these quantities in the porous medium 
and the dependent variables must be regarded as averages of the corresponding 
physical quantity over the cross-section of the capillary. 

The system of equations Eq. (|TJ) - (J2J) is closed by Poisson's equation of elec- 
trostatics, which in SI units is written as 



dE 



N 



e— = p e (x,t) = 2_^ez l c l (x,t). (3) 

i=i 

Here p e (x,t) is the charge density and e is the permittivity of the medium. Time 
variations in electrophoresis are sufficiently slow that any errors resulting from 
replacement of the complete Maxwell's equations of electromagnetism by Poisson's 
equation, Eq. ([3|), is negligible. Eq. (UJ) ((3j) provide a closed system of equations 
for calculation of the time evolution of the functions c, (x, t) if the initial conditions 
are known and if appropriate boundary conditions are specified. 



1 this relation is not valid for macro-ions or colloidal particles owing to Debye layer effects. 
The relation is valid provided the ionic size is much smaller than the Debye length 
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Fig. 1 A schematic sketch showing the CE setup and the electropherogram (inset) recorded 
by the detector. The diagram is not to scale. 



3 Electroneutrality and the Kohlrausch function 



If 'a' denotes a characteristic length scale over which the concentrations may be 
considered to vary (this may be taken, for example, as the peak width), then it may 
be shown that the ratio of the term on the left hand side of Eq. ([3j) to that on the 
right is of the order of (Ajj/a) , where Xp is the Debye length in the buffer. Since 
in CE, a is of the order of microns while Arj is in the nanometer range, (Arj/a) 
is extremely small. As a result, Eq. ([3]) may be replaced by the electroneutrality 
condition: 



pe{x,t) 



N 

E 



eZiCi(x,t) = 0, 



(4) 



first invoked by Planck [14] in the context of the problem of the liquid junction 
potential. This condition caused some confusion in the early history of the sub- 
ject, since it appeared that in the absence of the charge density, no electric field, 
or in our problem, perturbation to the applied field is possible. The controversy 
was resolved [19] when it was realized that Eq. ^ should be interpreted in an 
asymptotic rather than in an exact sense. 

The disappearance of the left hand side of Eq. ([3]) may at first sight appear to 
be problematic because we have now lost our equation for determining E. This is 
not in fact a difficulty as E can be determined directly from Ohm's law and the 
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condition of the constancy of current: 

J = J = aE (5) 

where the electrical conductivity iqj 

N N 

a = 22^i z i ec i = e 2 2J Zj^Cj (6) 

i=i i=i 

and Jo is a constant. The first of the two equalities in Eq ([5]) is a consequence of 
charge neutrality but the second requires certain additional assumptions. To see 
this, substitute g, from Eq. (J2j) in the expression for the current 



N d 

i=l 



N 



/ ZZiCiUi 



N 

+ Ey j mzjecj. (7) 



If the diffusive contributions are now neglected, then the second part of Eq. ([3]) , 
J = aE or Ohm's law follows. 

The constraints imposed by Eq. Q and ([5]) have the effect of replacing two 
of the differential equations governing the transport by algebraic relations. It was 
pointed out by Kohlrausch [2b that a further reduction of the governing equations 
is possible if the diffusion terms are neglected in the ion transport equations. To 
see this we define the Kohlrausch function: 

N 

K(x,t) = J2%- ( 8) 

t—1 

If we now take the derivative of both sides of this relation with respect to time 
and use Eq. (|TJ|, © and Q we get 



(9) 
i=x Ui ) 

If the right hand side is neglected we get 

K(x,t) = K(x,0), (10) 

and this provides an algebraic relation to replace one more differential equation 
for transport. 

In contrast to the electro-neutrality condition which is respected to high ac- 
curacy, the neglect of the diffusion terms in arriving at either Ohm's law or the 
Kohlrausch relation is questionable in applications related to CE. The ratio of the 
diffusive term to the other terms is measured by the inverse of the Peclet number, 
Pe = va/D in terms of the characteristic peak width a and velocity v. While it is 
true that the Peclet number is generally large in the applications we are concerned 
with here, the transit times L/v across the length (L) of the capillary is also very 
long. Thus, the ratio of the diffusive broadening (y / '2DL/v) to a characteristic 




2 the dependance of conductivity on concentrations is actually nonlinear at high concentra- 
tions 20 . However we ignore this effect as it does not fundamentally alter our analysis except 
in replacing one nonlinear dependence of ion velocity on concentration by another 
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Fig. 2 Sketch illustrating the basis of the assumption that the Kohlrausch function may be 
assumed constant in the sample zone shortly after sample injection. 



peak width (a) is measured by the square root of (L/a)Pe~ which is generally 
not small even though Pe ~ 10 — 100. Furthermore, when the profile steepens, 
as it does in electromigration dispersion, the second derivative term multiplying 
Di in the transport equations can become very large so that the diffusive term 
once again has a finite effect no matter how small the diffusivity is. Nevertheless, 
the Kohlrausch relations have been fruitfully applied for certain special situations 
such as in calculating fluxes across interfaces separating two homogeneous regions. 
In the next section we show how the ideas of Kohlrausch and Ohm's law can be 
generalized to treat the case of nonzero diffusivities. 



4 Reduction to a one equation model 



Let us now consider an idealized model of electrophoresis where only three species 
are involved: a sample or analyte ion, a positive ion and a negative ion. In place 
of the suffix i we will from now on use a suffix 'p' to denote the positive ion, 'n' 
to denote the negative ion and no suffix at all to indicate the sample. Thus, z p 
is the signed valence of the positive ion, z n that of the negative ion and z that 
of the sample. Clearly z p is positive, z n is negative and z could be of either sign 
(but not zero, since a neutral species does not electromigrate). We will further 
assume that all the species have the same mobility (w), and therefore the same 
diffusivity D = ukgT, though their electrophoretric mobilities m = eziU could 
be widely different. Many analytes of importance in electrophoresis are linear 
polyelectrolytes, and for these molecules the charge scales linearly with the mass 
M while the diffusivity scales roughly as M^ 1 ' 3 . It might therefore be argued 
that diffusivities of ions vary over a relatively small range when compared to 
their electrophoretic mobilities, though it must be acknowledged that our primary 
motivation for introducing this assumption is the pragmatic one of wanting to 
reduce the mathematical complexity of the problem. 
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It is readily seen that by virtue of our assumption Di = D and the elec- 
troneutrality condition Eq. Q, the first term of Eq. (0 vanishes so that Ohm's 
law, Eq ([S]) is recovered. Further, in view of our assumption of equal diffusivities, 
Eq. ([9]) reduces to the diffusion equation 



f=°e <»> 



which admits an analytical solution: 



K(x,t) = ^== / tf(y,0)exp 



-t-oo 

dy. (12) 



1 f ~ „. ~ (x-y) 



ADt 
If Koo denotes the unperturbed value of K, then the perturbation in K, SK = K 



Koo spreads diffusively from the injection region so that its width, w ~ y w o + 2Z)£, 
where wq is the width of the peak at injection. The essential observation that allows 
us to extend Kohlrausch's ideas to the diffusive case is illustrated in Figure [2] and 
it is the following: even though K(x,t) is in general a time dependent quantity, 
in the vicinity of the analyte peak we may assume it to be a constant equal to 
its unperturbed value K = Koo ■ This is because the distance between the analyte 
peak and the injection zone (X) increases linearly with time X ~ fiE^t where Eq is 
the applied field so that at long times, X ^> w. The minimum length of time that 
needs to pass for this to happen is max [wo/v, D/v 2 ] where v = fj,Eo is the velocity 
of an isolated analyte ion in the back ground electrolyte and max [ ] indicates 
the larger of the two numbers within the brackets. It is instructive to express this 
in terms of the diffusion time r^ = a 2 jD the time it takes the analyte to diffuse 
over a characteristic distance. If we identify this characteristic distance with wo, 
that is, a = wo, then the condition may be written as tjr^ ;§> max[Pe~ , Pe - ] 
where Pe = va/D is the Peclet number introduced earlier. If we use the estimates 
v ~ 1 mm/s, a = wq ~ 10 /im, D ~ 10~ 5 cm 2 /s which are fairly typical, we have 
Pe ~ 10 and r d ~ 0.1 s. Thus, the condition K = Koo in the vicinity of the analyte 
peak may be used with confidence for times t 2> 10 -2 s. This is indeed very small 
in comparison to the time interval between injection and detection which is tens 
of minutes. In simple terms, SK spreads so slowly that the analyte peak leaves it 
behind essentially in the time it takes to traverse a distance of the order of its own 
width. 

Setting the right hand side of Eq. ([8j) to the constant value Koo immediately 
yields a third algebraic relation in addition to Eq. (0} and Eq. ([5| between the set 
of four dependent variables: c p , c n , c and E. If we use these to express all of the 
remaining variables in terms of c and substitute in the transport equation for c, 
we obtain, after some algebra that we omit, a single nonlinear partial differential 
equation describing the evolution of the analyte concentration c. This equation 
may be conveniently expressed in terms of the dimensionless quantity <f> = c/c^ 
where c„ is the concentration of the negative ions in the background electrolyte: 

m +v di-{i—^) =D d^- (13) 

The dimensionless parameter a is defined as 

a = (z - Z n )(z - Zp) = (fl - n n ){n - ftp) 
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Since \ip > and /j, n < 0, a is positive if fi is in the interval {^ n ,^,p) and is 
negative otherwise. The parameter a is essentially a dimensionless form of the 
"velocity slope parameter" that has been used [22] to characterize electromigration 
dispersion, since av = V'(0) where V{4>) = v/(l — a<f) is the local value of the 
migrational velocity of an analyte ion within the sample zone. 

The concentrations of the negative and positive ions may be readily expressed 
in terms of 4> using the charge neutrality condition, Eq Q and the defining relation 
for K(x,t), Eq (JSJ : 



c n — Cn 






1+ " 



Zn 



z„ 



+ _0^SK, (15) 

Zp — Zfi 



M ™ -5K. (16) 



Zn 



Zp Zp Zn 

It is seen from the above equations that once the analyte has moved out of the 
injection zone (0 — > 0) the background ions do not immediately return to their 
equilibrium concentrations. A fluctuation is created in the concentration of the 
background ions which does not propagate, but decays in place on a diffusive time 
scale. If the capillary carries an electroosmotic flow and if the sensor is sensitive 
to the background ions (e.g. a UV absorbance detector) these fluctuations may be 
detected with the analyte peaks complicating the interpretation of the signal. The 
existence of these anomalous peaks are well known [23 ,24,25 and are variously 
referred to as "water peaks" , "system peaks" or "vacancy peaks" . 

When a is positive, the evolution equation for cf>, Eq (|13|) has a singularity at 
4> = l/a which needs to be understood. Clearly, this singularity arises because 
of the vanishing of the conductivity a in Ohm's law Eq ((5}. Since m = ezjit, 
the vanishing of the conductivity when <f> — > l/a implies that one or both of 
the background ion concentrations must have become negative. This, of course is 
physically impossible and would imply that at least one of our two assumptions, 
namely charge neutrality or constancy of the Kohlrausch function would have to be 
violated before this singularity is reached. In fact, the highest analyte concentration 
to which Eq (|13|) may still be used can be calculated by requiring that both 
the conditions c n > and c p > be respected where c n and c p are given by 
Eq. (|15[) and Eq. (|16[) with 8K = 0. These inequalities may be reduced after some 
algebra to the requirement <p < cj> c where <j> c = (z p — z n )/(z p — z) when z < and 
4> c = — [z n (zp— Zn)/[zp(z — Zn)] when z > 0. As a practical example, if z = ±2z p and 
z n = —z p we get 4> c = 2/3 which would correspond to conditions of extremely high 
sample loading. If initial conditions are used such that (f> > 4>c in some regions, 
there will be an initial transient period that is not described by Eq (|13[) , but once 
diffusion causes the analyte concentration to drop below the critical value, the 
subsequent evolution should follow Eq. (|13[) . 



5 Analysis of the one equation model 



Eq. (|13[) in general describes a one dimensional nonlinear wave [26J. The nonlin- 
earity arises because the local wave speed V{4>) = v/(l — a<fr) changes with 4>. In 
particular, if a > 0, larger values of <j> correspond to higher speeds, so that if <j> 
has a maximum at some point, then analyte ions near the maximum out runs the 
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slower ions ahead of it causing the distribution of <j> to lean forward and develop a 
sharp gradient at its leading edge (a shock wave). An analogous situation arises if 
a < except in this case the distribution of (f> leans backwards as seen in Panels C 
and D of Figure 3. The diffusive term on the right hand side of Eq. (JT3J) works to 
counter the sharpening of the wave due to the nonlinearity. When it is dominant 
the wave approaches a Gaussian shape in time with small departures from this 
limiting shape due to the nonlinearity. In case of small diffusivity the diffusive 
term is small everywhere except in the vicinity of the shock; its presence prevents 
the sharp change in 4> a t the shock from evolving into a true mathematical dis- 
continuity. In the remainder of this section we attempt to rationalize some of the 
observed facts about electromigration dispersion by exploiting some well known 
results from the mathematical theory of nonlinear waves. 



5.1 The weakly nonlinear case 

Let us suppose that \a<j>\ <C 1 for all x and t. Since </> is the ratio of analyte ion 
to that of the background ions, typically \<f>\ <C 1. The above condition can also 
be satisfied even when \4>\ is not small compared to unity, provided \a\ is small, 
that is, if the electrophoretic mobility of the analyte ion closely matches that 
of one of the background electrolytes. In either case, one may wish to neglect 
the term a<j> in the denominator of Eq. (|13(1 . which would then describe waves 
that propagate with constant velocity v while simultaneously undergoing diffusive 
spreading. This of course describes the usual "linear" electrophoresis. To see the 
effect of electromigration dispersion we must treat the term a<j> as small, but not 
zero. This is achieved by making the approximation (1— q</>) _1 « l+a<j> in Eq. (|13l) . 
which then becomes: 

d(j> . d(t> dcf> d 2 4> . _. 

Eq. (|17[) takes a more familiar form if we describe it in a reference frame translating 
with velocity v, that is, we transform variables from (x,t) to (£,i) where £ = x — vt. 
Eq. (|17p then becomes the one dimensional Burgers' equation^! 






^+2*v^ = D^. (18) 



Eq. (jl8(l is well known in fluid mechanics and is studied as a model for break- 
ing water waves, one dimensional shocks in gas dynamics and the like. It can be 
reduced [26] to the one dimensional diffusion equation through a nonlinear trans- 
formation of variables and thus solved exactly. The solution is: 

/•+oo p _ / r+oa 

*(£,*) = / |^|exp{-.P(»j; £,*)}*?/ / exp{-P(rnZ,t)}dn (19) 



3G 



where 



P fae.*) = %r^ + TT / <K£',0K'. (20) 



3 in contrast to Eq. 1131 ) - the solutions of Burgers' equation are well behaved for all <f>, though 
the solution would correspond to physical reality only if < <j> < <j>c 
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Given any initial profile for </> its subsequent time evolution may be calculated 
simply by evaluating the integrals in Eq. (|19[) . A special case of interest is 



<j>(x,0) =rs( x ) (21) 

where F = J_°° </>(x,t)dx (a conserved quantity) and 5(x) is the Dirac delta func- 
tion. This may be an appropriate model for t 2> t<j when the peak width is much 
larger than its initial value so that the detailed shape of the initial concentra- 
tion distribution may be disregarded. We will use the quantity r to characterize 
the degree of sample loading. It has dimensions of length and may be attributed 
the following simple meaning: it is the length along the capillary that contains 
the same number of negative ions as there are analyte ions in the injected peak. 
Substitution of Eq. ([21]) in Eq. ([19]) yields 

, 1/2 „ fx~vt" 



^=s(s) 'm) ™ 



where 

Ap-x 2 

F ^ = l + { A/2)er H x) <& 

with erfc(x) = (2/*/ir) J exp(— t ) dt denoting the complementary error function. 
The constant A is given by A = e — 1 where P = Fv/D is a Peclet number based 
on the characteristic length scale F. If aP is small, A « and in this case the de- 
nominator in Eq. (|23[) may be approximated by one. In this limit Eq. (|22[) reduces 
to a spreading Gaussian profile as expected. If aP is significant, A is non-negligible 
and it has the same sign as a. In this case the denominator in Eq. (|23|) creates 
an asymmetry in the profile, Eq. (|22[) . causing it to steepen in the direction of 
peak propagation or create a steep trailing edge depending on the sign of a. This 
behavior of course is well known from experimental and numerical studies of elec- 
tromigration dispersion and is often referred to as "fronting" or "tailing" . It may 
be of interest to note that Eq. (|23[) is similar to a certain empirical function that 
has historically been used to fit experimental data on asymmetric peaks [27U28] . 
However, unlike this so called "Haarhoff-Van der Linde (HVL) function" Eq. ()23[) 
is not an empirical fit but a mathematical consequence of the governing equations. 
The quantity \a\P that determines the strength of this effect may be expressed in 
terms of several relevant time scales: \a\P = r e T,i/T a where r a = a/v is an advective 
time, Td is the diffusive time scale introduced earlier and r e = \a\F/v is a new time 
scale associated with sample loading. The importance of the scale r e will become 
apparent later, but here we provide some numerical estimates. If we assume that 
the analyte peak is a Gaussian with a 10 /im width and a peak concentration that 
is 2 percent of the background, then r e ~ 0.5 ms, if v ~ 1 mm/s and \a\ ~ 1. 
The other two timescales using the estimates provided earlier are r a ~ 10 ms and 
r d ~ 100 ms. Thus, in this example, P ~ 0.5 and therefore A ~ 1 which would 
correspond to "high" sample loading and on the basis of Eq. (|23|) strong deviations 
from Gaussianity may be expected. Thus, the quantity \a\P may be regarded as 
a suitable measure of sample loading in the context of peak distortion. 

Another useful quantity that may be calculated [26] from the theory of non- 
linear waves is the time needed to develop a shock: 



2av | 

ox 



(24) 
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where the suffix 'min' indicates the minimum value of the function (which is clearly 
negative for a unimodal distribution) . If the initial profile is a Gaussian of standard 
deviation wo then the minimum is easily evaluated and we get 



*. 



/7re w 



o 



2 avr 






(25) 



If we substitute our numerical estimates, we get t s hock ~ 0.2 s which is very short 
compared to the total transit time across the capillary, typically tens of minutes. 
Shock like solution profiles may be obtained [26] from Eq. (|22|) on taking the 
asymptotic limit of D — > 0. If a > we have leading edge shocks: 



<*'*> = Ta 



x h 
x f 



vt. 



x b 



±-1 

x b 



1 + 2 



for xt > x > x b , 



(26) 
(27) 
(28) 



and if a < 0, we have trailing edge shocks: 



(j>(x,t) 

x f 
x b 



1 

2Jo[ 

vt, 
x f 



X 

x f 



for xt>x>x b , 



(29) 
(30) 
(31) 



Here Xf and x b are the locations of the "front" and "back" of the distribution 
4>{x,t) which vanishes outside the interval [aJi,,a;/]. 

The profile, Eq. (|22p . may be used to calculate some useful quantities such as 
the location of the centroid: 



t c {t) 



J_ x(f>(x,t) dx 
L™<i>(x,t)dx 



(32) 



and the variance: 



Kt) 



f_™(x-x c ) 2 <j)(x,t) dx 



J-oo <f>(.x,t) dx 



(33) 



The evaluation of the integrals is facilitated on transforming to the new variable 
r) = (x—vt)/v4JDt. The second moment about the centroid {x c ), Eq (|33[) . is related 
in a simple way to ctq which is defined as in Eq (|33[) but with x c replaced by xq = vt. 
This relation (sometimes referred to as the "parallel axis theorem" in mechanics) : 
op = (xq — x c ) 2 + u 2 easily follows from Eq (|33p . Its use simplifies the calculation 
of the second moment, and we get after some algebra: 



x c (t) = vt + V4Dt^-, 

-TO 



and 



a 2 (t) = ADt 



El-El 

F F 2 



(34) 
(35) 
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where P*. = /_ °° r\ F(rj) dr\ is the k th moment of the function P(r?) defined by 
Eq. (|23|) . In general the P& need to be determined numerically except for k = 0, in 
which case Fq = aP^/ii which is merely a statement of the fact that the solution, 
Eq (|22p . satisfies the normalization implied by Eq (|21|) . The speed of the centroid 
V c = dxc/dt readily follows from Eq (|34[) : 



V c {t) 



1± El .fll 



3 V £ 



(if P^oo). (36) 



with the + or — sign depending on whether a is positive or negative. If the effective 
diffusivity is defined by o\ (t) = 2-D e ff t then 



.-■2 



D eff = 2rW||-|U~±|a|P V (if P^oo). (37) 



The asymptotic forms for large P following the symbol ~ in Eq. (|36|1 and Eq. (|3T|) 
are easily obtained on calculating the moments Pj. using the profiles Eq. (|26|l and 
Eq. (|29[) . Thus, we see that the centroid velocity approaches v asymptotically for 
times that are large compared to r e , however the spread of the peak is diffusive at 
all times and may be described by the effective diffusivity indicated by Eq (|37p . 
The phrase "all times" here refers to times that are nevertheless large compared 
to r d since otherwise Eq (|22[) may not be applied and the solution in this regime 
depends on the detailed features of the initial profile. 



5.2 The strongly nonlinear case 

The assumption \a<f>\ <C 1 may not be a good one especially during the early part 
of the evolution in the case of high sample loading. The replacement of Eq. (|13|) 
by Burgers' equation, Eq. (|18[) . is not a good approximation in that case and 
we refer to this as the strongly nonlinear regime. However, an analytical solution 
is still possible if Pe 3> 1. In this case, the diffusion term in Eq. ()13|) may be 
neglected except in the immediate vicinity of a shock. Since Pe is usually large 
in electrophoresis it is reasonable to treat the problem as one of propagation of a 
shock front for times t > t s hock! which formally corresponds to setting the right 
hand side of Eq. (|13p to zero while admitting discontinuous functions within the 
class of solutions. It should be pointed out that Eq. (JT3J) , with D = 0, and its 
solution Eq. (|42[) had been noted earlier by others [6,8, though its consequences 
for electromigration dispersion had not been adequately worked out. We do so 
here. 

A solution to Eq. (|13[) with D = may be souglrqj in similarity form in terms 
of the variable X = x/vt, that is, we look for solutions of the type 4>{x, t) = 'P(X) 
where <f (X) is an unknown function. Substitution of this form in Eq. (|13[) with D 
set to zero results in the following requirement that <P(X) needs to satisfy: 



d 
dX 



(p 



1 — a<3> 



d$ 
X dX' ^ 



4 the existence of such similarity solutions is suggested by the invariance of the equations 
under a scale transformation: x — > x/p and t —} t/p for any nonzero constant p 
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which may be easily rearranged as 



X 



(l-a$y 



d$ 

dX 



0. 



(39) 



Thus, <P must either satisfy the differential equation d<P/dX = and therefore must 
be a constant or an algebraic equation that is easily solved: 



*(*) = I 



1± 



(40) 



It is clear that neither of these functions by itself can represent the solution for all 
X but rather each describes a piece of it. The branch <P is constant must describe 
the solution outside the analyte zone and because of the requirement $(±oo) = 
the constant must be zero. The solution within the analyte zone is best discussed 
by considering the cases corresponding to the two signs of a separately. 



5.2.1 Leading edge shocks (a > 0) 



The solution branch Eq. (|40[) can be matched at one of its boundaries to the branch 
$(X) = only if we accept the negative sign. This boundary then corresponds to 
X = 1 which must be identified with the trailing edge. The leading edge is at 
some unknown position X = X* > 1. We can determine Xt by exploiting the 
condition that the integral of cf> is a conserved quantity, which immediately follows 
on integrating Eq. (|13|) : 



r 



+ oo j-Xf 

<j){x,t) dx = vt I <P(X)dX 
-oo J\ 



It is best to present the solution in terms of our original variables: 



4>(x,t) = - 
a 

Xf, = vt. 
Xf = vt 



I _ /£* 



Wt 



for xr>x> x^, 



(41) 

(42) 
(43) 
(44) 



and 4>(x,t) = outside the interval [xf,,xr]. Here r e = \a\T/v is the time scale 
introduced earlier. Two quantities of particular importance in electrophoresis are 
the maximum value of ij>, 4>m{i) = cf>(xf,t) and the peak width w = xj — x\>. These 



(45) 
(46) 

(47) 



may be easily calculated from Eq. (|32 

(t>m(t) = - 



1 + 



U)(t) = Vs/Tet 

The propagation speed of the peak is 



2 ■ J 2. 



Vf 



d.v 



./■ 



dt 
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It is easy to verify by direct substitution that Eq. (|45[) and Eq. (|47|) satisfy the 
Rankine-Hugoniot condition for shock waves [26] . 

It is clear from these formulas that in the asymptotic limit of long times (( > t e ) 
the wave propagates unchanged in shape with its width increasing in proportion 
to the square root of time and its amplitude decreasing in inverse proportion to 
the square root of time while its integral is conserved. This is classical diffusive 
behavior, and therefore an effective diffusivity (-D e ff) can be introduced through 
the asymptotic relation a ~ 2D e fj t for the increase of variance with time. To 
determine D e ff one needs to calculate the variance of the distribution represented 
by Eq. (|42D and then determine the leading order term in its asymptotic expansion 
in the small parameter r e /t. This may be done and leads after some algebra to the 
rather simple final result 

D eS = ~\a\rv. (48) 

We write \a\ in place of a in anticipation of the fact that this expression is equally 
valid for a trailing edge shock (a < 0). This is seen to be identical to the large 
Peclet number (P) limit of the effective diffusivity obtained from the solution of 
Burgers' equation indicated in Eq. (137 



5.2.2 Trailing edge shocks (a < 0) 



The requirement that $ must be non-negative implies that only the negative sign 
in Eq. (|40p corresponds to physically realizable solutions. Thus, 



*(*) 



k.v| 



l 



- l 



(49) 



where X < 1. Therefore in this case, X = 1 actually corresponds to the leading 
edge and X = X b < 1 is the trailing edge which is where the shock is located. To 
determine X b we once again use the constancy of the integral of </>: 



r 



<f>(x,i) dx = vt / ${X)dX. 

•oo JX b 



The relations analogous to Eq. (f4"2)) - (|4"4")l are 



(50) 




for Xf > x > Xf,, 



(51) 



x b = vt I- \j -j- (52) 

Xf = vt. (53) 

The amplitude, width and speed of the peak may be calculated as before, they are 

, -l 



4>m{t) 



kvl 



)(t) = v^/r e t 



1 



t 



T/ dx b 



l- 



(54) 

(55) 
(56) 
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The effective diffusivity D e g is once again given by Eq. 

Since Burgers' equation follows from Eq. (|13p if (j> is small, one might expect 
that Eq. (|51[) - (|53p should reduce to Eq. (|29p - (|3ip at asymptotically large times. 
This is indeed true. To show this, first observe that from Eq. (|52|) 



Xf, = vt 



1 



vt 



vt 



1 



(57) 



at large times. Secondly, if we put x = x t — £ and regard (,/xt as small (since 
w ~ \ft while Xf ~ t), then 



4>(x,t) 



J_ 

M 



j_ 

M 



-1/2 



e 



2 1 CV I ;r 



/ 



1 

2U 



x 

x f 
(58) 



Thus, Eq. ([29]) - ([STJ are recovered. Similarly, if a > 0, Eq. (g2j) - (gj) go over to 
Eq. (PHI) - (J2BJ) at large times. 



6 Numerical simulations 



In order to test the accuracy of some of the results presented in the previous 
section, they were compared to solutions of Eq. (|13[) computed numerically using 
a second order (central) finite difference algorithm for spatial derivatives coupled 
with a fourth order Runge-Kutta time stepping scheme with fixed grid and step 
size. We will refer to this numerical calculation as "the full solution" . 

As a test problem, an initial profile (/>(x,0) = 4> m exp(— x 2 /2a 2 ,) was specified 
and 4>m was varied keeping <jq fixed. This involves no loss of generality because if 
Eq. (| 13p is rewritten in terms of dimensionless space and time variables: x* = x/ao 
and £* = vt/ao, then it is easily seen that the only parameters that enter into the 
problem are: a, P = vF/D and F/ao = V2TV(j> m - Comparisons were made over a 
wide range of parameter values, but for brevity, we present here only a few cases 
that are illustrative of the general trend. 

In Figure [3] we compare the analyte concentration profiles obtained by nu- 
merically evolving Eq (|13|) with the analytical formulas Eq. (|22[) and Eq. ([42j) or 
Eq. (|5ip . It is seen that for large values of P all of the profiles are very close. For 
smaller value of P, the full solution is in good agreement with the analytical solu- 
tion of the Burgers' equation, Eq. (|22p . but does not agree with the shock profiles 
Eq. (|42p or Eq. (|5ip . This is because the dynamics is dominated by diffusion. 

The formulas given by Eq. (|37p and Eq. (|48p could be of some practical value 
so we compared them with the corresponding results from the full solution in 
Figure|U To do this, the variance a 2 was calculated from the concentration profile 
in the full solution and the rate of change of variance d{a 2 )/dt was plotted as 
a function of time t. In all cases, this quantity was found to asymptote to a 
constant value, and equating this value to 2D e g- gave an effective value of the 
diffusivity. These are shown in Figure [3] by the symbols for a range of values of the 
Peclet number P, for a = ±0.5. It should be observed that P e g is independent of 
the sign of a as predicted by the theoretical results. The small difference in the 
value of D e R for positive and negative a is due to the fact that the peak width 
approaches its asymptotic value from above or from below depending on the sign 
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Fig. 3 Comparison of the normalized concentration <f> (y-axis) as a function of (x — vt)/(To 
(x-axis) obtained from the full solution, Eq. (1 13 I t (open circles), Eq. H22II (solid lines) and 
Eq. ||42]| or Eq. f5TT) (dashed lines) at fixed instants of time vt/ao = 50, 100 and 200. The 
parameters are: P = 62.7: a = 0.5 (Panel A), a = -0.5 (Panel C) and P = 5.0: a = 0.5 (Panel 
B), a = -0.5 (Panel D). 



of a as is suggested by Eq. (|46p and Eq. (|55l) for the peak width. Thus, the 
numerical determination based on the full solution results in a slight over estimate 
of the asymptotic value in one case and a slight under estimate in the other. The 
numerical calculation is seen to be in good agreement with Eq. (|37|) . However, 
the dashed line has a constant offset with respect to the solid line because when 
P is large, D e ff/D = 1 + 5M-P + ■•■, and the limit represented by Eq. (|48[) is 
obtained by neglecting the first term in relation to the second. It is therefore a 
valid approximation and has a low relative error when P is much larger than unity, 
but if this is not the case then one must use Eq. (I3T 



7 Concluding remarks 



In this paper a simple model is constructed and rigorously solved with the objective 
of providing a physical understanding of the essential features observed in the 
CE signal at high sample loading. The model leaves out many features that are 
undoubtedly present in any real experiment, such as, the presence of many more 
than three ionic species, shifts in ionic equilibria, differential diffusion between 
species and the like. On the other hand, the fact that the principal qualitative 
features of the signal may be reproduced in this highly simplified model indicates 
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Fig. 4 The normalized effective diffusivity as a function of the Peclet number P. The solid line 
is from Eq. H37H and the dashed line is from Eq. H48II . The symbols indicate values calculated 
from the numerical solution of Eq. I13I I. Two cases are shown, a = +0.5 (circles) and a = —0.5 
(triangles). 



that the neglected effects, though present, are not causative but only incidental to 
the problem. 

The main results that follow from this analysis that might have some validity 
in a broader context are: 

1. Leading edge shocks are to be expected when the analyte ion has an elec- 
trophoretic mobility that is intermediate to that of the co and counter ions 
and a trailing edge shock should be seen otherwise. 

2. A peak showing a leading edge shock will elute earlier and a peak showing a 
trailing edge shock will elute later compared to the arrival times in an otherwise 
identical experiment where sample loading is kept low. 

3. The quantity jajP, where P = Pv/D is a Peclet number based on the length 
scale r, completely determines the nature of the wave. If \a\P << 1 the peak 
would evolve into a Gaussian with a small amount of asymmetry and if \a\P >> 
1 a shock like structure is formed. It is the single controlling parameter that 
must be kept low in order to mitigate the effects of electromigration dispersion. 

4. At long times the peak spreads diffusively with an effective diffusivity D e g. The 
ratio D e fi/D depends on the Peclet number P; it is proportional to P when P 
is large. 

5. The peak asymmetry that is characterized by the quantity A has an exponential 
dependence on the Peclet number, P: A = e — 1, as a result, the transition 
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from a "weak asymmetry" to a "strong shock" as one increases the sample 
loading is very sharp. 
6. If the analyte electrophoretic mobility is closely matched to the electrophoretic 
mobility of one of the background species (this makes \a\ small) peak distortion 
and dispersion can be minimized. 

These observations appear to be consistent with laboratory experiments 14 , 29 and 
results from numerical simulations [7J. The analytical formulas presented in this 
paper can be checked quantitatively through carefully designed experiments that 
closely mimic the simplifying assumptions made in this study. It is hoped that this 
paper will motivate efforts in that direction. 
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